This document was prepared for a briefing to Dr. Angell on analyzing the equity implications of Public Safety Power Shutoffs (PSPS) in California.

In Lieu of data on the actual locations and times of Californians affected by PSPS - which are expecting from PG&E soon, we examined the overlay of CalFIRE’s Fire Hazard Severity Zones (FHSZ) with two (2) statewide datasets describing health equity:

This preliminary analysis, with its limitations, highlights a few key findings:

  1. Inclusive health equity indices may not best describe disparate vulnerabilities in areas likely affected by PSPS.
  2. The racial makeup of the highest risk areas is predominantly white.
  3. Far more adults with a disability live in the highest risk areas than in the lower risk areas.
  4. High fire risk near urban areas tends to be among wealthier individuals than in rural and frontier Medical Services Study Areas.

Fire Hazard Severity

The CalFIRE data is a decent proxy for areas potentially affected by PSPS. The state identifies Moderate, High, and Very High risk areas of both state and local responsibility (can be explored in map below). There are other datasets of wildfire risk that we could explore that may include areas of federal responsibility.

pal <- colorFactor(palette = "Reds",levels = c("Moderate","High","VeryHigh"),ordered = T,
                        domain = fhsz_sp$FRSZ,
                        n = 3)


pal2 <- colorQuantile(palette = c("#2172B4","#9DCAE1","#A4DBA1","#008836"), n=4,  domain = 0:100)


maphpi <- merge(tracts_sp, { 
          
          hpi[,.(ct10, hpi2_pctile, fire_risk, city, county_name)]%>% mutate(city = ifelse(city =="","[No City Listed]",city))
  
          
        })



hpimap2 <- filter(maphpi, hpi2_pctile <= 25 & fire_risk == "VeryHigh") %>% mutate(city = ifelse(city =="","[No City Listed]",city))

# here we create the Leaft map with two layers, one for all ages and one for elderly
  leaflet() %>%
     # addProviderTiles(providers$Esri.WorldStreetMap, group = "Street Map") %>%
     #  addProviderTiles(providers$Stamen.Terrain, group = "Terrain") %>%
     #  addProviderTiles(providers$Esri.WorldImagery, group = "Sattelite") %>%
    addProviderTiles(providers$Stamen.Toner, group = "B/W") %>%
    addPolygons(data = fhsz_sp,
        color = "#444444",
        weight = 1,
        smoothFactor = 0.1,
        fillOpacity = 0.8,
        fillColor = ~ pal(fhsz_sp$FRSZ),
        stroke = FALSE,
        label = ~fhsz_sp$FRSZ,group = "Fire Hazard Severity Zones"
      ) %>%  
        addPolygons(data = maphpi,
          
        color = "#444444",
        weight = 1,
        smoothFactor = 0.1,
        fillOpacity = 0.8,
        fillColor = ~ pal2(maphpi$hpi2_pctile),
        stroke = FALSE,
        popup = paste0("This location is in the top ",round(maphpi$hpi2_pctile,0),"the percentile of all tracts in the state for the healthy places index"),
        group = "HPI"
      ) %>% 
        addPolygons(data = hpimap2,
        color = "#444444",
        weight = 4,
        smoothFactor = 0.1,
        stroke = TRUE,
        popup = paste0("This tract has 'Very High Fire Risk' and is among the most vulnerable (bottom 25%), according to HPI. It is near ", hpimap2$city,", ", hpimap2$county_name," County"),
        group = "HPI + Fire Risk") %>% 
    addLayersControl(
    # baseGroups = c("Street Map", "Terrain", "Sattelite", "B/W"),
      overlayGroups = c("Fire Hazard Severity Zones", "HPI", "HPI + Fire Risk"),
      
                    options = layersControlOptions(collapsed = FALSE)
                ) %>%
      setView(lat = 37.085206, lng = -119.540085, zoom = 6) %>% #defaults the view of the original map to show the entire state
      addEasyButton(easyButton(
        icon="fa-globe", title="Zoom to State",
        onClick=JS("function(btn, map){ map.setView([37.085206, -119.540085],6); }"))) %>% # adds a button to return to the whole state view
      addLegend("bottomleft",
                pal = pal,
                values = factor(c("Moderate","High","VeryHigh"), levels = c("Moderate","High","VeryHigh")),
                opacity = 1,
                title = "Fire Hazard Severity Zones"
                ##lables = c("X%-X% (most vulnerable)","X%-X%","X%-X%","X%-X%","X%-X% (least"%. In this tract, the "mapTemp$def," is higher than "mapTemp),))##
      ) %>%  # adds a button to return to the whole state view
      addLegend("bottomleft",
                pal = pal2,
                values = na.omit(maphpi$hpi2_pctile), 
                labels = c("bottom 25%","","","healthiest 25%"),
                opacity = 1,
                title = "HPI Percentile") %>%
                ##lables = c("X%-X% (most vulnerable)","X%-X%","X%-X%","X%-X%","X%-X% (least"%. In this tract, the "mapTemp$def," is higher than "mapTemp),))##
      
    leaflet::hideGroup("HPI")
  DT::datatable({

  
    foo <-maphpi %>% select(county_name, city, ct10, hpi2_pctile, fire_risk) %>%
      rename(County = county_name,
             City = city,
             CensusTractID = ct10,
             HPI_percentile = hpi2_pctile, 
             Fire_Risk = fire_risk)
    
    st_geometry(foo) <- NULL
    
    foo
    
    },  filter = 'top',
    options=list(pageLength = 5)
  ) %>% DT::formatRound(4, 0) %>% DT::formatString(3)

Health Equity Index and Fire Risk

By combining the Fire Hazard zones with the health equity data, we can begin to see how populations living in these different areas compare. For example, using the HPI, we can compare the average Index score (and components) for tracts in the different Fire Risk Severity Zones. The actual scores are not very meaningful, but we can see that we see increase HPI Scores for tracts in the highest FHSZs.

foo <- hpi[,.(ct10, fire_risk, hpi2score, economic, education, housing, healthcareaccess, neighborhood, pollution, transportation, social)] %>% 
  melt.data.table(id.vars = c("ct10","fire_risk"), variable.name = "hpi_component",value.name = "value")

hcboxplot(var = foo$hpi_component, x = foo$value, var2 = foo$fire_risk,
          outliers = FALSE) %>% 
  hc_chart(type = "column") %>% 
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "HPI Index and Components")) %>%
  hc_yAxis(title = list(text = "HPI Score (higher = healthier)")) %>%
  hc_title(text = "Healthy Places Index and Fire Hazard Severity Zones",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from CalFIRE and Public Health Alliance of Southern CA",
              align = "left", style = list(fontWeight = "bold"))

Elderly (by race) and Fire Risk

A couple of the indicators in the CCHVIs may be more illustrative of the concentration of vulnerability in the highest fire risk areas. For example, the number of elderly individuals (over 65) that live in Very High Risk FHSZs is greater than the number in Moderate and High FHSZs combined, and the majority are white.

goo <- cchvi_tracts[ind == "elderly" & race !="Total", .(number = sum(numratr, na.rm = T)), by = .(fire_risk, strata, race)] %>%
  .[,race := factor(race, levels = c("AfricanAm","AIAN","Asian","Latino","White","Multiple","Total"))] %>% na.omit()


hchart(goo, "bar",  hcaes(group = fire_risk, y = number, x = race)) %>%
  hc_add_theme(hc_theme_smpl())  %>% 
  hc_xAxis(title = list(text = "Adults over age 65")) %>%
  hc_yAxis(title = list(text = "Race")) %>%
  hc_title(text = "Elderly by Race and Fire Hazard Severity Zones",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from CalFIRE and CDPH",
              align = "left", style = list(fontWeight = "bold"))

Poverty (by race) and Fire Risk

Examining the number of individuals living in poverty (<200%FPL) by race in Very High Risk FHSZs highlights more of the non-white vulnerability in these areas.

goo <- cchvi_tracts[ind == "poverty" & race !="Total", .(number = sum(numratr, na.rm = T)), by = .(fire_risk, strata, race)] %>%
  .[,race := factor(race, levels = c("AfricanAm","AIAN","Asian","Latino","White","Multiple","Total"))] %>% na.omit()


hchart(goo, "bar",  hcaes(group = fire_risk, y = number, x = race)) %>%
  hc_add_theme(hc_theme_smpl())  %>% 
  hc_xAxis(title = list(text = "Individuals in Poverty (<200% FPL)")) %>%
  hc_yAxis(title = list(text = "Race")) %>%
  hc_title(text = "Poverty by Race and Fire Hazard Severity Zones",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from CalFIRE and CDPH",
              align = "left", style = list(fontWeight = "bold"))

Disability and Fire Risk

Another indicator in the CCHVIs of interest is individuals with a disability. The number of individuals living with a disability (mental or physical) is greatest in the highest risk fire areas.

joo <- cchvi_tracts[ind == "disability" , .(number = sum(numratr, na.rm = T)), by = .(fire_risk, strata)] 


hchart(joo, "column",  hcaes(group = fire_risk, y = number, x = strata)) %>%
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "Type of Disability")) %>% 
  hc_yAxis(title = list(text = "Number of Adults with Disability")) %>%
  hc_title(text = "Number of Adults with Disability in Fire Hazard Severity Zones",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from CalFIRE and CDPH",
              align = "left", style = list(fontWeight = "bold"))

Percentage of People in Poverty

Income affects the ability of households to deal and cope with many of the ‘inconveniences’ of climate change. We see the mean percentage of people living in poverty to be very different between Very High fire risk areas in the urban vs. rural and frontier areas of the state (using Medical Services Study Area definitions).

goo <- cchvi_tracts[ind == "poverty" & race =="Total"]


hcboxplot(var = goo$MSSA_DEFINITION, x = goo$est, var2 = goo$fire_risk,
          outliers = FALSE) %>% 
  hc_chart(type = "column") %>% 
  hc_add_theme(hc_theme_smpl()) %>% 
  hc_xAxis(title = list(text = "Urbanicity - Medical Services Study Area")) %>%
  hc_yAxis(title = list(text = "Percent in Poverty (<200% FPL)"), min = 0, max = 100) %>%
  hc_title(text = "Tract Level Poverty by Fire Hazard Severity Zone and Urbanicity",
           margin = 20, align = "left") %>%
  hc_subtitle(text = "data from CalFIRE and CDPH",
              align = "left", style = list(fontWeight = "bold"))